Dataset

I use the fire incident dispach from opendata of NewYork

Loading Packages

library(dplyr)
## Warning: il pacchetto 'dplyr' è stato creato con R versione 4.2.3
## 
## Caricamento pacchetto: 'dplyr'
## I seguenti oggetti sono mascherati da 'package:stats':
## 
##     filter, lag
## I seguenti oggetti sono mascherati da 'package:base':
## 
##     intersect, setdiff, setequal, union
library(ggplot2)
## Warning: il pacchetto 'ggplot2' è stato creato con R versione 4.2.3
library(tidyverse)
## Warning: il pacchetto 'tidyverse' è stato creato con R versione 4.2.3
## Warning: il pacchetto 'tibble' è stato creato con R versione 4.2.3
## Warning: il pacchetto 'tidyr' è stato creato con R versione 4.2.3
## Warning: il pacchetto 'readr' è stato creato con R versione 4.2.3
## Warning: il pacchetto 'purrr' è stato creato con R versione 4.2.3
## Warning: il pacchetto 'stringr' è stato creato con R versione 4.2.3
## Warning: il pacchetto 'forcats' è stato creato con R versione 4.2.3
## Warning: il pacchetto 'lubridate' è stato creato con R versione 4.2.3
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ lubridate 1.9.3     ✔ tibble    3.2.1
## ✔ purrr     1.0.2     ✔ tidyr     1.3.0
## ✔ readr     2.1.4
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(lubridate)
library(mapview)
## Warning: il pacchetto 'mapview' è stato creato con R versione 4.2.3
library(sf)
## Warning: il pacchetto 'sf' è stato creato con R versione 4.2.3
## Linking to GEOS 3.9.3, GDAL 3.5.2, PROJ 8.2.1; sf_use_s2() is TRUE
library(geojsonio)
## Warning: il pacchetto 'geojsonio' è stato creato con R versione 4.2.3
## Registered S3 method overwritten by 'geojsonsf':
##   method        from   
##   print.geojson geojson
## 
## Caricamento pacchetto: 'geojsonio'
## 
## Il seguente oggetto è mascherato da 'package:base':
## 
##     pretty
library(leaflet) 
## Warning: il pacchetto 'leaflet' è stato creato con R versione 4.2.3
library(broom)
## Warning: il pacchetto 'broom' è stato creato con R versione 4.2.3
fire_data <- read.csv("Fire_Incident_Dispatch_Data_last_50k.csv")
alarm_box_loc <- read.csv("Alarm_Box_Locations.csv")
firefighter_stations <- read.csv("fdny-firehouse-listing.csv")
fire_data$ALARM_BOX_NUMBER <- as.character(fire_data$ALARM_BOX_NUMBER)
fire_data$INCIDENT_BOROUGH <- as.factor(fire_data$INCIDENT_BOROUGH)
 
fire_data <- fire_data %>% rename(al_number = ALARM_BOX_NUMBER, borough = INCIDENT_BOROUGH)

firefighter_stations$Borough <- as.factor(firefighter_stations$Borough)
firefighter_stations <- firefighter_stations %>% rename(borough = Borough)

alarm_box_loc_new <- alarm_box_loc %>% select(BOROBOX, LATITUDE, LONGITUDE) %>% 
                        rename(al_number = BOROBOX, lat = LATITUDE, lon = LONGITUDE)

fire_data <- fire_data %>%
  mutate(al_number = str_pad(al_number, width = 4, pad = "0"))
fire_data_new <- fire_data %>%
  mutate(al_number = case_when(
    borough == "BRONX" ~ paste0("X", al_number),
    borough == "BROOKLYN" ~ paste0("B", al_number),
    borough == "MANHATTAN" ~ paste0("M", al_number),
    borough == "QUEENS" ~ paste0("Q", al_number),
    borough == "RICHMOND / STATEN ISLAND" ~ paste0("R", al_number),
    TRUE ~ al_number
  ))
fire_data_merged <- merge(fire_data_new, alarm_box_loc_new, by = "al_number", all.x = TRUE)

fire_data_merged <- na.omit(fire_data_merged)
firefighter_stations <- na.omit(firefighter_stations)
alarm_inc_count <- fire_data_merged %>%
                    group_by(al_number, lat, lon) %>%
                    summarise(Count = n()) %>%
                    rename(incident_count = Count)
## `summarise()` has grouped output by 'al_number', 'lat'. You can override using
## the `.groups` argument.
alarm_inc_count_sdf <- st_as_sf(alarm_inc_count, coords = c("lon", "lat"), crs = 4326)
count_inc_brough <- count(fire_data_merged, borough)
count_inc_brough <- count_inc_brough %>% rename(incident_count = n)
count_inc_brough <- count_inc_brough %>% mutate(borough = recode_factor(borough, "BRONX" = "Bronx", "BROOKLYN" = "Brooklyn", "MANHATTAN" = "Manhattan", "QUEENS" = "Queens", "RICHMOND / STATEN ISLAND" = "Staten Island"))
count_inc_brough
##         borough incident_count
## 1         Bronx           4882
## 2      Brooklyn           7867
## 3     Manhattan           7094
## 4        Queens           5756
## 5 Staten Island           1350
firefighter_stations <- firefighter_stations %>% rename(lat = Latitude, lon = Longitude)
firefighter_stations_sdf <- st_as_sf(firefighter_stations, coords = c("lon", "lat"), crs = 4326)
head(firefighter_stations_sdf)
## Simple feature collection with 6 features and 10 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -74.01252 ymin: 40.70347 xmax: -73.9929 ymax: 40.71976
## Geodetic CRS:  WGS 84
##                                              FacilityName       FacilityAddress
## 1                                      Engine 4/Ladder 15       42 South Street
## 2                                     Engine 10/Ladder 10    124 Liberty Street
## 3                                                Engine 6     49 Beekman Street
## 4 Engine 7/Ladder 1/Battalion 1/Manhattan Borough Command  100-104 Duane Street
## 5                                                Ladder 8 14 North Moore Street
## 6                                       Engine 9/Ladder 6       75 Canal Street
##     borough Postcode Community.Board Community.Council Census.Tract     BIN
## 1 Manhattan    10005               1                 1            7 1000867
## 2 Manhattan    10006               1                 1           13 1075700
## 3 Manhattan    10038               1                 1         1501 1001287
## 4 Manhattan    10007               1                 1           33 1001647
## 5 Manhattan    10013               1                 1           33 1002150
## 6 Manhattan    10002               3                 1           16 1003898
##          BBL
## 1 1000350001
## 2 1000520022
## 3 1000930030
## 4 1001500025
## 5 1001890035
## 6 1003000030
##                                                                           NTA
## 1 Battery Park City-Lower Manhattan                                          
## 2 Battery Park City-Lower Manhattan                                          
## 3 Battery Park City-Lower Manhattan                                          
## 4 SoHo-TriBeCa-Civic Center-Little Italy                                     
## 5 SoHo-TriBeCa-Civic Center-Little Italy                                     
## 6 Chinatown                                                                  
##                     geometry
## 1 POINT (-74.00754 40.70347)
## 2 POINT (-74.01252 40.71007)
## 3 POINT (-74.00525 40.71005)
## 4 POINT (-74.00594 40.71546)
## 5 POINT (-74.00668 40.71976)
## 6  POINT (-73.9929 40.71521)

Downloand of the GEO JSON

geojson_newyork <- geojson_read("NYC_BoroughBoundaries.geojson",  what = "sp")
geojson_newyork <- setNames(geojson_newyork, c("borough_code", "borough", "shape_area", "shape_leng"))
geojson_newyork$borough <- as.factor(geojson_newyork$borough)
geojson_newyork$borough_code <- NULL
head(geojson_newyork)
## class       : SpatialPolygonsDataFrame 
## features    : 5 
## extent      : -74.25559, -73.70001, 40.49613, 40.91553  (xmin, xmax, ymin, ymax)
## crs         : +proj=longlat +datum=WGS84 +no_defs 
## variables   : 3
## names       :       borough,    shape_area,    shape_leng 
## min values  :         Bronx,  1187174772.5,  325917.35395 
## max values  : Staten Island, 636520502.758, 888199.731385
geojson_newyork@data = data.frame(geojson_newyork@data, count_inc_brough[match(geojson_newyork@data$borough, count_inc_brough$borough),])
geojson_newyork@data$borough.1 <- NULL
mapview(list(firefighter_stations_sdf, geojson_newyork),
        zcol = list(NULL, "incident_count"),
        legend = list(FALSE, TRUE),
        homebutton = list(FALSE, TRUE), layer.name = list(NULL, "indicents_number"), alpha.regions = 0.5, aplha = 1)
mapview(list(alarm_inc_count_sdf, geojson_newyork),
        zcol = list("incident_count", "borough"),
        legend = list(TRUE, TRUE),
        homebutton = list(TRUE, TRUE),
        layer.name = list("incident_count", "brorough"),
        alpha.regions = list(0.7, 0.4), 
        aplha = list(1, 1))